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Abstract 

In this paper, we study a flower population in which self-reproduction is not permitted. In- 
dividuals are diploid, that is, each cell contains two sets of chromosomes, and distylous, that is, 
two alleles, A and a, can be found at the considered locus S. Pollen and ovules of flowers with the 
same genotype at locus S cannot mate. This prevents the pollen of a given flower to fecundate its 
own stigmata. Only genotypes AA and Aa can be maintained in the population, so that the latter 
can be described by a random walk in the positive quadrant whose components are the number of 
individuals of each genotype. This random walk is not homogeneous and its transitions depend on 
the location of the process. We are interested in the computation of the extinction probabilities, as 
extinction happens when one of the axis is reached by the process. These extinction probabilities, 
which depend on the initial condition, satisfy a doubly-indexed recurrence equation that cannot be 
solved directly. Our contribution is twofold : on the one hand, we obtain an explicit, though intri- 
cate, solution through the study of the PDE solved by the associated generating function. On the 
other hand, we provide numerical results comparing stochastic and deterministic approximations 
of the extinction probabilities. 

Keywords: Inhomogeneous random walk on the positive quadrant; boundary absorption; transport 
equation; method of characteristics; self-incompatibility in flower populations; extinction in diploid 
population with sexual reproduction 

AMS: 60G50; 60J80; 35Q92; 92D25 

1 Introduction 

We consider the model of flower population without pollen limitation introduced in Billiard and Tran 
[4]. The flower reproduction is sexual: plants produce pollen that may fecundate the stigmata of other 
plants. We are interested in self-incompatible reproduction, where an individual can reproduce only 
with compatible partners. In particular, self-incompatible reproduction prevents the fecundation of 
a plant's stigmata by its own pollen. Each plant is diploid and characterized by the two alleles that 
it carries at the locus S, which decide on the possible types of partners with whom the plant may 
reproduce (as it encodes the recognition proteins present on the pollen and stigmata of the plant). 
We consider the distyle case with only two possible types for the alleles, A or a. The plants thus 
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have genotypes AA, Aa or aa. The only interesting case is when A is dominant over a (see [4]), and 
we restrict to this case in this work. Then, the phenotype, i.e. the type of proteins carried by the 
pollen and stigmata, of individuals with genotypes AA (resp. Aa and aa) is A (resp. A and a). Only 
pollen and stigmata with different proteins can give viable seeds, i.e. pollen of a plant of phenotype A 
can only fecundate stigmata of a plant of phenotype a and vice- versa. It can be seen that seeds AA 
cannot be created, since the genotype of individuals of phenotype a is necessarily aa that combine only 
with individuals of phenotype A that have genotypes AA or Aa, therefore we can consider without 
restriction populations consisting only of individuals of genotypes Aa and aa. Each viable seed is 
then necessarily of genotype Aa or aa with probability 1/2. It is assumed that ovules are produced 
in continuous time at rate r > and that each ovule is fecundated to give a seed, provided there 
exists compatible pollen in the population. The lifetime of each individual follows an exponential 
distribution with mean 1/d, where d > 0. In all the article, we consider 



r > d 



(1.1) 



which, we will see, is the interesting case. 
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Figure 1: (a) Transition rates for the continuous-time pure-jump Markov process (Xt,Yt)teM.+ - 
(b) Transition probabilities of the embedded random walk, that we denote by (Xt,Yt)t<=m (here and 
throughout, N is the set {0, 1,2,.. .} ), with an abuse of notation. 



Let us denote by X t and Y t the number of individuals of genotype Aa (phenotype A) and aa 
(phenotype a) at time t £ M+. The process (Xt, Yt) t £R + is a pure-jump Markov process with transitions 
represented in Fig. 1(a). A stochastic differential equation (SDE) representation of (X t , Y t ) te ^ + is given 
in [4]. Here we forget the continuous-time process, and we are interested in the embedded discrete- 
time Markov chain, which we denote, with an abuse of notation, by (Xt, Yj)t g N, and with transitions 
represented in Fig. 1(b): 

PyPi^i) = (i- i,i)] = - ( - ,1' v ViA(*i,Yi) = (* + 

(r + d)(i+j) 

Pijp^n) = (i,j - 1)] = {r+ d d) J {i+j y = (i,3 + 1)] 

where Pjj means that the process starts with the initial condition (Xq,Yq) = (i,j). The main and 
profound difficulty is that this random walk is not homogeneous in space, while techniques developed 
in the literature for random walks on positive quadrants mostly focus on the homogeneous case (see 
e.g. Fayolle et al. [6], Klein Haneveld and Pittenger [8], Kurkova and Raschel [9], Walraevens, van 
Leeuwaarden and Boxma [13]). We introduce a generating function (1.4) that satisfies here a partial 
differential equation (PDE) of a new type that we solve. Although the particularity of the problem 
is exploited, these techniques and the links between probability and PDEs may be extended to carry 
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out general studies of inhomogeneous random walks in cones. The introduction of PDEs through 
generating functions had been already used by Feller [7] for a trunking problem with an inhomogeneous 
random walk in dimension 1. To our knowledge, the case of inhomogeneous random walks in the cone 
with absorbing boundaries has been left open. In [10], the discriminatory processor-sharing queue is 
considered but boundaries are not absorbing and the overall arrival rate is constant, which is not the 
case in our model. 

When one of the phenotype A or a disappears, reproduction becomes impossible and the extinction 
of the system occurs. We are interested in the probability of extinction of (X t , Yt)teN (or, equivalently, 
in that of (Xt,Yt)teBL+)- Let us introduce the first time at which one of the two types gets extinct: 

t = ia£{t G N : X t = or Y t = 0}. (1.2) 

For i, j G N, let us denote by 



the absorption probabilities, and by 



Pi,j = Pijfo < oo] (1.3) 



P(x,y) = ^ Pijx'y* (1.4) 



their generating function. By symmetry arguments, we have, for all i,j G N, 

Pi,j=Pj,i- 0-5) 
Moreover, for any i, j G N such that i = or j = 0, we have 

Pij = 1. (1-6) 

In Section 2, we will see that the Pij's satisfy the Dirichlet problem associated with the following 
doubly-indexed recurrence equation 

di dj r r 

(r + d)(i + j) (r + d)(i + j) 2(r + a) 2(r + a) 

and with the boundary condition (1.6). This problem does not admit simple solutions. There is no 
uniqueness of solutions to this problem. Note that the constant sequence equal to 1 is a solution. 
However, we are interested in solutions that tend to as i or j tends to infinity, since, [4] (see Proposi- 
tion 2.2 in this paper), estimates for p^j were obtained through probabilistic coupling techniques; they 
show that in the case (1.1) we consider, is strictly less than 1. In fact, the Pij's correspond to the 
smallest positive solution of the Dirichlet problem, and are completely determined if we give the prob- 
abilities (pi i)i3si- We conclude the section with more precise estimates of the absorption probabilities 
Pij as the initial state (i, j) goes to infinity along one axis (Proposition 2.3). These new estimates rely 
on Proposition 2.2 and on comparisons with one-dimensional random walks. In Section 3, we consider 
the generating function P(x,y) associated with the Pij's and show that it satisfies a PDE, that has 
one and only one solution, that is computed (Proposition 3.5) explicitly with a dependence on the 
(pil)i^i, prompting us to use the name "Green's function". This provides a new formulation of the 
solution of (1.7), that is however uneasy to work with numerically. Hence, in Section 4, we propose 
two different approaches leading to numerical approximations of the solution of the Dirichlet problem 
(1.6)— (1.7), that are based on stochastic and deterministic approaches. 

In conclusion, we provide here several approaches to handle the extinction probabilities of the in- 
homogeneous random walk (Xt,Yt)t£R + of our problem. Estimates from [4] are recalled and the 
recurrence equation (1.7) is solved numerically and theoretically, pending further investigation of the 
PDE formulation. 
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2 Existence of a solution 
2.1 Dirichlet problem 

We first establish that the extinction probabilities Pij's (1.3) solve the Dirichlet problem (1.6)-(1.7). 

Proposition 2.1. (i) The extinction probabilities (pij)i,i>l are solutions to the Dirichlet problem 
(1.7) with boundary condition (1.6). Uniqueness of the solution may not hold, but the extinction 
probabilities (pij)ij£N define the smallest positive solution to this problem. 

(ii) Let the probabilities (p» 1)^1 be given. Then the probabilities (Pi,j)i,j^i are completely determined. 



Proof. We begin with Point (i). Equation (1.7) is obtained by using the strong Markov property at 
the time of the first event. Let us denote by K the transition kernel of the discrete-time Markov chain 
{X t ,Y t )tm*; we have: 



Kf(iJ) 



(/(* + 1,3) + f(iJ + 1))^7-K + f(i ~ U) , , t. , ^ + f(i,3 ~ 1)- 



2(r + d) JK ,J ' (r + d){i+j) JK,J '{r + d){i+j)' 

Following classical proofs {e.g. [3, 11]), the extinction probabilities (pij)tjeN satisfy the equation: 

Vi.jeN*, f(i,j) = Kf(i,j) and V*,j € N, /(*, 0) = f(0,j) = 1. (2.1) 

The constant solution equal to 1 is a solution to (2.1). Let us prove that (pij)ij£^ is the smallest pos- 
itive solution to (2.1). Let / be another positive solution. Let us consider M t = /p4 n f{t,-n)}> ^nf{t,T })j 
with To defined in (1.2). Denoting by (Gt)teN the filtration of (Mt)teN, we have: 

E[M t+ i | Gt] = E[M t+1 l T0<t + M t+l l T0>t | Gt] 

= E[M t t T0<t + f(X t+l ,Y t+1 )l T0>t | Gt] 
= M t l T0<t + t T0>t E[f{X t+1 ,Y t+1 ) | Gt] 
= M t l T0<t + l T0>t Kf(X t ,Y t ) 
= M t l T0<t + l T0>t f(X t , Y t ) = M t . 

Hence (Mate's is a martingale, which converges on {to < oo} to f(X TQ ,Y T0 ) = 1 (see the boundary 
condition in (2.1)). Thus by using the positivity of / and Fatou's lemma, we obtain that for every 
i,j G N: 

f(i,j)=Eij[Mo] = EmEy[M t ] ^ E[limmf M t l T0<oo ] = E itj [l n<00 ] = PiJ . 

This concludes the proof of Point (i). 

Let us now consider Point (ii). Assume that the probabilities (pi,i)i>i are given, and let us prove, 
by recursion, that every pij can be computed. By symmetry, we only need to prove that this is the 
case for i ^ j. Assume 



(Hrec j): for j € N* all the pk,e J s for I ^ j and k ^ I can be computed from the Pi i's 
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and let us prove that we can determine the pij+i's for i ^ j + 1. From (1.7) we get: 

2(r + d) 2di 2d j 

ftj+i = ^ PiJ ~ r(»+j) Pf " 1J ~ Ki + i)^' -1 ~ J ' 

All the terms in the r.h.s. of (2.2) are known by (Hrec j), and hence Pij+i can be computed for any 
i ^ j + 1. This concludes the recursion. ■ 

The following result shows that there is almost sure extinction in the case r ^ d. In the interesting 
case r > d, it also shows that there is a nontrivial solution to the Dirichlet problem (1.6)-(1.7). 



Proposition 2.2 (Proposition 9 of [4]). We have the following regimes given the parameters r and d: 

(i) If r ^ d, we have almost sure extinction of the population. 

(ii) If r > d{> 0), then there is a strictly positive survival probability. Denoting by (i,j) the initial 
condition, we have: 




In Point (ii), only bounds, and no explicit formula, are available for the extinction probability p^j. 
The purpose of this article is to address (1.7) by considering the Green's function P{x,y) introduced 
in (1.4). 

2.2 Asymptotic behavior of the absorption probability as the initial state goes to 
infinity along one axis 

In this part, using the result of Proposition 2.2, we provide more precise estimates of the asymptotic 
behavior of the absorption probability pij = pj t i when j — > oo. In particular, these estimates will be 
very useful when we tackle the deterministic numerical simulations (see Section 4.2). 



Proposition 2.3. If j — > oo, then 

2dl 2d(r 2 + dr + 2d 2 ) 1 / 1 \ 

PIJ = PjA = — " 21 _^A\ ^2 + h3 ' ( 2 - 3 ) 

r j r z (r + a) j z \3 J 

Proof. In addition to to, defined in (1.2), we introduce 

5 = inf{t G N : Y t = 0}, T = wf{t € N : X t = 0}, 

the hitting times of the horizontal axis and vertical axis, respectively. Note that we have to = inflS 1 , T}. 
Let / : N — > N be a function such that /(j) < j for any j ^ 1. In the sequel, we will choose f(j) = [ej\ , 
with e G (0, 1) and where [-J denotes the integer part). We obviously have the identity: 

pij = P (lj) [r < oo] = P (1)i) [r < f(j)] + r (1J) [f(j) < T < oo]. (2.4) 

To prove Proposition 2.3, we shall give estimates for both terms in the r.h.s. of (2.4). 

First step: Study of PnjJro ^ f(j)}- Since /(j) < j, it is impossible, starting from (1, j)> to reach the 
horizontal axis before time /(j), and we have P^j^to ^ /(?)] = ^(ij)^ ^ /(?')]■ or der to compute 
the latter probability, we introduce two one-dimensional random walks on N, namely X~ and X + , 
which are killed at 0, and which have the jumps 

F i [Xf = i-l]=qt, F l [X±=i + l]=pt, W> i [xt=i]=rf ) q f+pf+rf = l, 



where 

± _ di ± ^ r 

Qi (r + d)(i + j^f(j)Y Pi 2(r + d)' 1 • ' 

Both X~ and X + are (inhomogeneous) birth-and-death processes on N. These random walks are 
implicitely parameterized by j. If = inf{i £ N : Xf 1 = 0}, then 

< /(j)] < P (1 j)[r < f(j)] < Pi[r + < /(j)]- (2.6) 

The quantities P^T 1 * 1 /(j)] are computable: we shall prove that 

Pl[T ± < f( 7 )l = ^ 1 2^ + ^ + 2^) 1 / 1 \ 

U " /UjJ r(i T /(i)) r*(r + d) (j qp f(j)) 2 \(j T /(?))/ ' 1 J 

The main idea for proving (2.7) is that the ^ being very small as j -4 oo, the only paths which will 
significantly contribute to the probability PifT 1 * 1 ^ f(j)] are the ones with very few jumps to the right. 
Let us define 

A^(p) = {the chain makes exactly p jumps to the right between and t} 
= {there exist 0^qi<---<q p ^t — 1 such that 

X fl + 1 ~ X qi = " ' = X q P + l ~ X q P = 1 ^ 

We are entitled to write 



PiP* «S f(j)} = Pi[^ < fU), A±. } (0)] +Pi[T± < f(j), A± y) (l)] 

+ Pi[r ± ^/(j),U p5j2 A± (i) (p)], (2.8) 

and we now separately analyze the three terms in the right-hand side of (2.8). First: 

fU) m 

F 1 [T ± < fU), A% } (0)] =^P 1 [T ± = k, A± w (0)] = £ (rf) fc ~V 

fc=i fc=i 

=r^(l-(rt) m ) =T^±(l + 0((rt) fU) )). (2.9) 
1 — r 1 — r i 

A Taylor expansion of qf /(I — r 1 ) according to the powers of q= f(j)) together with the fact that 
( r ± )/(i) = (l/(j =p /(j)) 3 ) provides that: 

2d ( 1 + ¥ „ / 1 



We now consider the second term in the right-hand side of (2.8). On the event A^r,.s(l), X^ first 
stays a time k\ at 1, then jumps to 2, where it remains fc 2 unit of times; it next goes to 1, and, after 
a time ks, jumps to 0. Further, since ^ /(j), we have fci + 1 + &2 + 1 + k$ + l ^ /(j)- Denoting 
by ki = k\ + /C3 the time spent in position 1, we thus have: 

Pipr^/OO.A^a)] 

E (»f)V(»t)* a ^ : (rf)*"^ : =pf^ E (^(^ 

fcl+fc 2 +fc3</(i)-3 fci+fc2</0')- 3 

P ^ 2± -(l + 0((rf Vr±) /(i) - 2 )). (2.11) 



(l-rf)(l-r* 



As for the first term, using the fact that (r^ V r^Y^ 2 = o(l/(j =p /(j)) 3 ) and a Taylor expansion 
according to the powers of 1 / ( j =F f(j)) gives that: 

ftP* < /«>, A^l)] . r(r + d ;f T/(j))2 (l + O (^ m )) . (2,2) 

Finally, let us consider the third term Pi[T ± ^ f(j), ^p>2^frj\(p)]- On L) p ^2-^f^(p), the two first 
jumps to the right are either from 1 to 2 and 2 to 3, or twice from 1 to 2. Thus, extinction means that 
there is at least 3 jumps from 3 to 2, 2 to 1 and 1 to or two jumps from 2 to 1 and 1 to 0. Since qf 
is an increasing function of i, we deduce: 

F 1 [T ± ^f(j),U p>2 Af U) (p)] 

From (2.8), (2.9), (2.11) and (2.13), we obtain (2.7). 
Second step: Study of P(ij)[/(j) < to < oo]. 

P( 1 ,,)[/(i)<T0<Oo] 

= E F adfd) <*> < °o\{x m ,Y fU) ) = (k,e))F {1J) [(x m ,Y m ) = (k,£)\ 

k,£^l 

= E p^ (1>j) [(x m ,Y m ) = (*>*)], ( 2 - 14 ) 

k,£^l 

by using the strong Markov property. Introduce now a function g : N — > N such that f(j) + g(j) < j 
for any j 1. We can split (2.14) into 

E Pu*M[(XfV)>YfV)) = M]+ E p^ { i^ x m^ Y m) = ^)\- ( 2 - 15 ) 

k,t>g(j) g{j)>k^l 

and / or 

With Proposition 2.2 we obtain the following upper bound for the first sum in (2.15): 

k,£^g(j) 

In particular, if we choose g such that as j —> 00, g(j) —> 00 fast enough, then clearly the term above 
is negligible compared to (2.7). For the second sum in (2.15), 

E PM p (U)[( J f/W. y /y)) = (*^]= E P^ m [(X m ,Y m ) = (k,£)}, 
g(j)^k>l g(j)>k^l 

and / or 

g(j)>m 

since by assumption j — f(j) > g(j) so that Y cannot reach values i ^ g(j) in f(j) steps. Then we 
have 

<V(M)[0^X m ^g(j)]. (2.16) 
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To obtain an upper bound for (2.16) we are going to use, again, a one-dimensional random walk. 
Introduce X, a random walk on N which is killed at 0, homogeneous on N* with jumps 

F k [X 1 = k-l]=q, F k [Xx = k + l]=p, F k [X 1 = k]=r, q + p + r = l, 

where 



* (r + d)(l + j-2f(j)y P 2(r + dY 
This walk is again parameterized by j. By construction of X, we have 

P (1J) [0 ^ X f{j) ^ g(j)} W x [X m < g(j)}. 

Denoting by m and a 2 the mean and the variance of (X2 — Xi), respectively (they could easily be 
computed), we can write 



X f{j) ~ m f^) < 9(j) - 1 - mf{j) 



By a suitable choice of the functions / and g, for instance f(j) = [ej\ with e G (0, 1) and g(j) = [j 3 ^\ , 
the central limit theorem gives that the latter is negligible compared to Fn ,j)i T o ^ fU)]- F° r this last 
term, using (2.7), (2.6) and letting e tend to provides (2.3). The proof is concluded. ■ 

Let us make some remarks on possible extensions of Proposition 2.3. 

Remark 1. 1. The proof of Proposition 2.3 can easily be extended to the asymptotic of pij as j — > 00, 
for any fixed value of i. In particular, we have the following asymptotic behavior: 

2. It is possible to generalize (2.8) by 

fc-i 

PiP* < SH)] = 5>i[ T± *S f(j), + Pi[T ± ^ /(j), U^A^]. (2.18) 

p=0 

We can show as in the proof of Proposition 2.3 that Pi < /(j), L> p>k A f p {j) } < (^ +1 ) fc+1 = o(l// +1 ) 

(see (2.13)) and that the probabilities F^T* < /(j), A^ (j) ] admit Taylor expansions in powers of 
T /(i)) wri ere the development for the p th probability has a main term in l/(j =p f(j)) p+1 - This 
allows us to push the developments in (2.7) to higher orders. 

For instance, the next term in (2.3) can be obtained by a long computation. First, we generalize 
(2.11) by writing that 

Pi[T± ^ f(j), 4 {j) ] 

= q^((pf) 2 (q^ 2 +pfpfqiq^) £ (rf)~ k > (rf)^ (2.19) 

fci+fc 2 +fc3^/(j)-5 

where k\, £?2 and k^ are the times spent by the random walk in the states 1, 2 and 3. Then, pushing 
further the Taylor expansion leads to: 

2d 1 2d(r 2 + dr + 2d 2 ) 1 

77-1 ■ = — 

rjTf(j) r 2 (r + d) (j =F f(j)) 2 

/2 r 2dv 2 24d(§+d) 5r 2 d 2 \ 1 / 1 \ 

+ Vr^ + r j r2(r + d) + 2(r + d) 2 f r +(i) 3j (jT/(j)) 3+ ^ {j ^ m y) ■ 



3 Green's function 



3.1 A functional equation for the Green's function 

In this section, we consider the Green's function P(x,y) defined in (1.4) associated with a solution 
of (1.7) in the same spirit as what can be found in Feller [7, Ch. XVII]. We show that it satisfies a 
non-classical linear PDE that can be solved (see Proposition 3.5). 



Proposition 3.1. (i) The function P(x,y) satisfies formally: 

AP(x,y) = h(x,y,P), (3.1) 



where: 



and where: 



dP dP 

AP(x, y) =Q(x, y)-^(x, y) + Q(y, x)^-{x, v) + R ( x > y)> ( 3 - 2 ) 



7" If* tjC 

Q(x,y) =(r + d)x - - - dx 2 , 

2 2 y 

r v 

R(yX " y ' = 2x + 2y ~ dX ~ dV " 



r 



"<■*• p) = - 2 <*• °> + rf ) + d *» ( ib + rb) ' (33) 

(ii) For given (pn)i^i, we have a unique classical solution to (3.1)-(3.3) on ]0, l[x]0, 1[. 



The function h in (3.3) only depends on a boundary condition (d 2 P/dxdy at the boundaries x = 
or y = 0, i.e. the p$ i's for i £ N*), which is non-classical, while the operator ^4 is of first order and 
hence associated with some transport equations. 

Proof of Proposition 3.1. Let us first establish (i). Using the Markov property at time t = 1: 

r . dj di 

2(r + d) (r + d)(i+j) (r + d)(i + j) 

then multiplying by x l y J , and summing over i, j G N* leads to: 
(r + d) ^2 (i + j)pi,jX % y 3 = - ^ (i + + PiJ+O^V 



2 



+ d J] •//'<•./ i + d ^ />; i ./•'' .'/' • (3.4) 



The l.h.s. of (3.4) equals: 

(d d \ 

x—+y—jP(x,y). (3.5) 
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For the r.h.s. of (3.4): 



j>0 j>1 



^l^PijC*- 1 ) 35 * V • .j ^ I'ij.i- 1 '' V 

i>2 



+ d^(j + l)ftjxV +1 + 1 ;/>,,.,•' V- 

For the first term in the r.h.s. of (3.6): 

J>1 J>1 J>1 J>1 

i>2 i>l i>l i>l 
r a rP(x,y) 

= 2Tx P ^ y) -2—^- (3 - 7) 

Similar computation holds for the 4th term of the r.h.s. of (3.6). For the 2nd term: 



j^2 j>2 i^l 

=i!| p ^ y) -^^'- (3 - 8) 

We handle the 3rd term of the r.h.s. of (3.6) similarly. Now for the 5th term: 

dJ2^ + 1)*V +1 =d(y 2 ^-P(x,y) + yP(x,y) + J>'y ) ■ ( 3 - 9 ) 

Similar computation holds for the last term of (3.6). From (3.5), (3.6), (3.7), (3.8) and (3.9) we deduce 
that: 

(r + <0 (4- + y) = ^ #-P(x, y) - + ^#P(x, y ) r P(X ' y) 



dx dy J 2 dx 2 x 2 dy 2 y 

ryd-., . r v-^ . ■ r x d _ , . r v-^ . 
+ 2x^ P(x ' y) " 2 + 2,^ /ir -' yi " 2^<"-'- 

+d^/ 2 ^P(x,y) + yP(x,y) + ^x l yj + d ^x 2 ^-P(x, y) + xP(x, y) + ^ , 
and finally J2i>iPiA ix ' 1 = x ^£^( x ^)- 

For point (ii), local existence and uniqueness stem from classical theorems [14]. Note that we construct 
an explicit, albeit complicated solution (3.5) using the method of characteristics. This concludes the 
proof. ■ 
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Remark 2. In the case of homogeneous random walks, the operator AP(x,y) has the product form 
R(x,y)P(x,y), see [9]. In some sense, this means that the inhomogeneity leads to partial derivatives 
in the functional equation. 

Remark 3. This technique allowing to compute the solution of a discrete, linear problem thanks to 
generating series is also called the Z-transform method. 

3.2 Characteristic curves 

In Sections 3.2 and 3.3 we establish, by the methods of characteristic equations, an explicit formula 
for the solutions to (3.1), which proves the existence of the solution. Since A is a first-order differential 
operator, we have a transport-like PDE. We introduce the following characteristic ordinary differential 
equations (ODEs). Let (x s ,y s ) s& ^ + be the solution to the system: 

x s = ^(s) = Q(x s ,y s ), 

dy (3 ' 10) 

Vs = -7-{s) = Q(y s ,x s ), 
v ds 

where Q has been defined in (3.2). The dynamical system (3.10) and its solutions will turn out to 
be decisive in the sequel — e.g. in Proposition 3.4, where we will use these characteristic equations in 
order to express the solutions to the fundamental functional equation (3.1). 



Proposition 3.2. For any initial condition (xo,yo) £ ^ 2 such that Xo,yo ^ 0, there exists a unique 
solution to (3.10), defined for s 6l by: 



with: 



Xre rs + [id e ds 
~ d(Xe rs + fie ds + 1)' 

2dx yo - d(x + y 



ds 



Xre rs + \ide 
d(Xe rs + ne ds - 1) ' 

2dx y Q + r(x + y ) 




^ S / ^ / / i I J 

s s / / / / / i 1-4 
//si s i t l J J 

////// 1 1 i ] 



i 1 1 i i 
i i 1. 1 1 ( j t i 
1 1 1 1 1 1 j i i 



(3.11) 



(3.12) 



Figure 2: Differential system (3.10): Vector field (Q(x,y),Q(y,x)). The plain lines correspond to 
the sets {Q(x, y) = 0} or {Q(y, x) = 0}. 

Proof. First of all, since Q(x, y) and Q(y, x) are locally Lipschitz-continuous for i,j/^0, there is local 
existence and uniqueness. Let us introduce the new variable z s = x s /y s , which is defined as long as 
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Figure 3: Differential system (3.10): (a) Neighborhood of the saddle point (1,1). (b) Neighborhood of 
the attractive equilibrium (r/d,r/d). 



y s / 0. Using the expression (3.2) of Q, the system (3.10) becomes: 

x s =(r + d)x s - -(1 + z s ) - dx 2 s , (3.13) 

. %sVs x s y s 

%s — n 
VI 

r z s r z s r z s r z s 

= (r + d)z s - dx s z s - (r + d)z s + h h dy s z s 

2x s 2y s 2y s 2x s 

=dx s (l-z s ). (3.14) 
From (3.14), we obtain: 



s d(i- Zs y 

and differentiating (3.14) with respect to the time s gives: 

dz s (l - z s ) + dz 2 



d 2 {l-z 
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Plugging this expression and (3.15) into (3.13) provides: 



dz s (l - z s ) + dz 2 s _ (r + d)z s r ^ + z ^ _ d % 



(3.15) 



d 2 (l-z s ) 2 d(l-z s ) 2 K d 2 (l-z s ) 2 ' 

and, therefore, 

z.(l - z s ) + 2z 2 - (r + d)(l - z s )z s + y (1 + z a )(l - z s ) 2 = 0. (3.16) 

The system (3.15)-(3.16) can be solved explicitly. Let us define u s = (1 + z s )/(l — z s ), so that 

u s - 1 



U s + 1 

from which we obtain: 

2u s .. 2u s (l + u s ) - 4u 2 s 



(3.17) 



/ii i2 ' ' - - i , r-i ' (3.18) 

(l + « s ) 2 (l + ii s ) d 
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Using these expressions in (3.16) provides: 



2u s (l + u a ) - Au 2 2 +2 _ Au 2 s 



(l + u s ) 3 u s + l (l + n s ) 4 

2 2u s rd 2u s 4 

-(r + d) + —- ^ =0 

v ; U S + 1(1 + U S ) 2 2« S + 1(!1 S |1) 2 

<=> 4« s (l + u s ) - 4(r + d)u s (l + u s ) +4rdu s (u s + 1) = 

«4> u s - (r + d)u s + rdu s = 0. (3.19) 



Hence it satisfies a second-order ordinary differential equation, that solves in: 

u s = Xe rs + fie ds , X,u€R. (3.20) 

From (3.15), (3.17) and (3.20): 

_ Xe rs + ue ds - 1 ^ Xre rs + ude ds 

Zs ~ Xe rs + ue ds + 1' Xs ~~ d(Ae rs + //e ds + 1) ' ^ ' ^ 

The integration constants A and u can be expressed in terms of the initial conditions xq and zq = xo/yo: 

_ 2dx -d(l + z ) _ r - 2dx + z r 

~ (l-z )(r-d) ' /i " (l- Z0 )(r-dY [6 - 22) 
This yields the announced result with y s = x s /z s . M 

Remark 4. Explicit expressions for x s and y s as functions of time and parameterized by r and d can 
be obtained. Using the expressions (3.12) of A and u in (3.21) and the relations between x s , y s and 
z s provides: 



xo + Vo- 2x y 

exp(rs) + exp(as) 



x + yo -2(d/r)x yo 



d x + y -2x yo , , , , , s , (1 - d/r)(y - x ) 

exp(rs) + exp(as) + 



(3.23) 



rx + y - 2(d/r)x y x + y - 2(d/r)x yo 

xo + yo- 2x y . . ... 

exp(rs) + exp(as) 



= x + j/o - 2(d/r)x y , g ^ 

d x + yo-2x y / \ . , , N (1 - d/r)(y - x ) 

exp(rs) + exp(as) 



rx + y - 2(d/r)x y x + y - 2(d/r)x y 

For the sequel, it is useful to notice that the constant 

A r xo + yo- 2x y 



(3.25) 



ud x + yo-2{d/r)x yo 

which appears in (3.23) and (3.24) belongs to ]0, 1[ for all (xo,yo) G]0, 1[ 2 ; this is due to the fact that 
d/r < 1. 

Below, we list some properties of the solutions to the dynamical system (3.10). The trajectories of 
the solutions to (3.10) can be decomposed into four steps. In order to describe them, let us introduce 



log 



Xo 



+ yo- 2x yo \ 



1 , ( V- <A \xo + yo-2(d/r)x y J 

S0 = — d log {- 17 J = d^r (3 - 26) 

and s± as the only positive roots of the denominators of x s and y s in (3.11): 

Xe rs + ne ds ± 1. (3.27) 
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Proposition 3.3. Let (xo,yo) G]0, 1[ 2 



(i 

(" 
(iii 

(iv 

The 

(v 
(vi 

(vii 
(viii 



The stationary solutions to (3.10) are the saddle point (1, 1) and the attractive point (r/d,r/d). 
so 6]0, oo[ and s± G]sq, oof. 

inf{s + , s_} = s + (resp. S-) if and only if yo < xq (resp. yo > xq). 

lim s ^ s x s = — oo and lim s ^ s _ y s = — oo. 

next points specify the behavior of the solution to (3.10). 

On [0, sq[, (x s ,y s ) belongs to ]0, 1[ 2 and goes to (x so ,y so ) = (0,0) as s — > sq. 

On ]s ,inf{s + ,s_}[, (x s ,y s ) goes decreasingly to (x- m ^ s+jS _y,y- m ^ s+iS _y)—by "decreasingly" we 
mean that both coordinates decrease. 

On ]inf{s + ,s_},sup{s+,s_}[, (x s ,y s ) goes to (x sup{s+ ^ s _ } ,y snp{s+rS _ } ) decreasingly. 
On ] sup{s + , s_}, oo[, (x s ,y s ) goes decreasingly to (r/d, r/d). 



Proof of Item (i). To find the stationary solutions to (3.10), let us solve x = and y = 0. With (3.10), 
we get Q(x,y) = and Q(y,x) = 0. This directly implies that (x,y) = (1,1) or (x,y) = (r/d, r/d), 
see Table 2. Let us now study the stability of these two equilibria. 
At (1, 1), the Jacobian of (3.10) is: 

v 

Jac(l, 1) = - J — dl, 

where J is the 2x2 matrix full of ones and where I is the 2x2 identity matrix. The eigenvalues of 
Jac(l, 1) are — d < and r — d > 0, associated with the eigenvectors (1, —1) and (1, 1), respectively. By 
classical linearization methods (e.g. [12, Chap. 3]), we deduce that the point (1, 1) is a saddle point. 
At (r/d,r/d), the Jacobian of (3.10) is: 

Jac(r/d, r/d) = — J — rl. 

The eigenvalues are — r < and d — r < 0, associated with the eigenvectors (1,-1) and (1,1), 
respectively. The point (r/d, r/d) is therefore attractive. ■ 

Proof of Item (ii) to (iv). Now we prove the different facts dealing with so, s+ and s_. First, (3.25) 
and the fact that r > d immediately imply that so E]0,oo[. Next, we show that (3.27) has on [0, oof 
only one root, which belongs to ]so,oo[. For this, we shall start with proving that (3.27) is positive on 
[0, so]- Then, we shall show that (3.27) is decreasing in ]so,oo[ and goes to — oo as s — > oo. 

In order to prove the first point above, it is enough to show that (3.27) is positive at s = and 
increasing on [0,so[. (3.27) is positive at s = simply because 

, , ... x + y 2y 2x (1 - d/r)[x + y ± (y - x )] 

A + n ± 1 = ± 1 = or = - > 0. 

yo -x Q y - xo yo - x x + yo~ 2(d/r)x yo 

To check that (3.27) is increasing on [0, so[, we note that the derivative of (3.27) is positive on [0, so[- 
actually by construction of so- 

Now we prove the second point. From (3.25) and since r > d, we obtain that (3.27) goes to — oo 
as s — > oo. Also, by definition of so, the derivative of (3.27) is negative on ]so,oo[, (3.27) is therefore 
decreasing on ]so, oof. 
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x 

Figure 4: Solutions corresponding to several initial conditions. We see that depending on the initial condition, 
the solutions converge of diverge to infinity. 



The fact that inf{s+, s_} equals s + (resp. s_) if and only if yo < xq (resp. yo > Xq) follows directly 
from (3.27). 

Finally, since the numerators of x s and y s are negative on ]sq,oo[, hence in particular at s±, it is 
immediate that lim.,^ x s = — oo and lim^ s _ y s = — oo. ■ 

Proof of Item (v) to (viii). Let us first consider Item (v). By definition of sq, the numerators of x s 
and y s in (3.23) and (3.24) vanish for the first time at sq. Moreover, since s± > so, both denominators 
are non-zero at so and x SQ = y SQ = 0. In particular, on [0,So[, we have (x s ,y s ) G In fact, 

(x s ,y s ) G]0, 1[ 2 . Indeed, on the segment {l}x]0, 1[, x s < whereas on ]0, l[x{l}, y s < 0: it is 
therefore not possible to go through these segments. 

We turn to the proof of Item (vi). Thanks to (3.23) and (3.24), just after the time so, (x s ,y s ) 
belongs to the negative quadrant M 2 . But for any (x, y) G ~S?_, Q(x, y) ^ and Q(y, x) ^ 0, see Table 
2, in such a way that both x s and y s are decreasing as soon as they stay in this quarter plane, in other 
words for s G]so, inf{s + , At time inf{s + ,s_}, one (or even the two if s + = s_, i.e. if xq = yo) 

of x s and y s becomes infinite. In the sequel, let us assume that inf{s+, s_} = s+; a similar reasoning 
would hold for the symmetrical case inf{s+, s_} = s_. 

Let us show Item (vii). Just after s+, (x s ,y s ) € (R + x M_)n{(x,y) G M 2 : Q(x,y) < 0, Q(y,x) < 
0}. The latter set is simply connected and bounded by the curve {(x,y) G M 2 : Q(x,y) = 0}, see 
Table 2. Using classical arguments (see e.g. [12]), we obtain that it is not possible to go through this 
limiting curve on which x s = 0; this is why for any s G]s+, s_[, (x s , y s ) remains inside of this set. 

We conclude with the proof of Item (viii). Just after the time s_, (x s ,y s ) G D {(x,y) G M 2 : 
Q(x,y) < 0, Q(y,x) < 0}. For the same reasons as above, (x s ,y s ) cannot leave this set and actually 
converges to (r/d, r/d). ■ 
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3.3 Use of the characteristic curves to simplify the functional equation 

Let us assume the existence of a solution P(x,y) to (3.1), and let us define g s = P(x s ,y s ). Then: 



dg . . dP dx s dP . dy s 

5s = ^ (s) = ^ (Xs ' ys) d7 + W (Xs ' ys) d7 



dP DP 

-Q^{x s ,y s )Q{x Sl y s ) + — (x s ,y s )Q(y s ,x s ). 



Thus, if P is a solution to (3.1), then: 

g s + R(x s ,y s )g s = h(x s ,y s ,P), (3.28) 

which looks like a first-order ODE for g, except that h depends on the boundary condition of P. 

We first freeze the dependence on the solution in h, i.e. we solve the ODE (3.28) as if the term 
in the right-hand side were a known function. Using the solutions to the characteristic equations, we 
shall prove the following result: 

Proposition 3.4. Let h(x,y) be an analytical function on [0,1[ 2 . Let (xo,yo) £ M 2 . The solution to 
the ODE 

g s + R(x s ,y s )g s = h(x s ,y s ), go = P(x ,y ), (3.29) 
where (x s ,y s ) s ^o are the solutions (3.11) to the characteristic curve starting at (xo,yo), is given by 

g h s =P(x , y ) exp (- j' R{x u , y u )du^ + j' h(x u , yu )e~ £ «(*-*«)d«d« 

=:F(s,x ,y ,h). (3.30) 

Proof. Equation (3.29) is an inhomogeneous first-order ODE. The solution to the associated homoge- 
neous equation is: 

g s = R(x ,y )exp f-J R(x u ,y u )du^j . 

The announced result is deduced from the variation of constant method. ■ 
A solution P to (3.1) hence satisfies the following functional equation for all s, xq and yo- 

P(x s ,y s ) =P(x ,y )e-foR^,vu)du + j* ^ ^ p)e -/J fl^^da^ (3 31) 

Jo 

with the function h defined in (3.3). Plugging the definitions (1.4) and (3.3) in (3.31), we obtain: 

m 

+d I" x uVu f — L_ + e- £ «(*-.Va)d« du . (3.32) 

Jo V 1 — x « >■ — y u J 

Notice that the r.h.s. of (3.32) depends only on the Pi i's and p± j's, while the l.h.s. depends on all Pi/s. 



r 
~2 
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Proposition 3.5. Let so > be defined in (3.26). We have: 

P(x ,y ) = JV K ,i fi^re/o^^^dn 

- d f x u y u ( — L_ + —LJ) e I u R{x a ,y a )da du _ (3-33) 

JO V 1 - x u 1-VuJ 

Before proving Proposition 3.5, let us show that the different quantities that appear in its state- 
ment are well defined — indeed, this is a priori not clear: as a — > sq, x a — > and y a — > 0, in such a 
way that R(x a ,y a ) — > oo, see (3.2). 

Lemma 3.6. Let i,j £ N. Then lim u ^ so (x u ) % (y u )i e-fo R(x a ,ya)da j s yj n ^ e if an d on \y ifi+j ^ 1 — and 
equals zero if and only if i + j ^ 2 . 



Proof. First, since the only zero of any function of the form aexp(as) + /3exp(6s) with a/3 < and 
a ^ b has order one, the following function has a simple zero at so : 

* / " n/ ' exp(ru) + exp(cki). (3.34) 



xo + yo- 2(d/r)x y 



Thanks to this and since s+, s_ > so, both x s and y s have a zero of order 1 at so- 
Moreover, with A and jj, defined in (3.12), we obtain that: 

X/fi + 1 - 



exp ( f R(x a ,y a )da , , , - 

VJo / (Vm) exp(ru) + exp(du) - l//i 

x — — — — ^— /dX/fx + exp((r + d)u). (3.35) 

(A///) exp(ru) + exp(du) + 1//J, (r / d){X/ n) exp(ru) + exp(du) 

It is indeed easy to check that the derivative of the logarithm of (3.35) is equal to R(x u ,y u ), for which 
we have an explicit expression, see (3.2), (3.23) and (3.24). 

From (3.35), we see that eh R( x a>ya)da three poles, namely at so,s+,s_. The zero at so has 
order one by using again the considerations on the zeros of (3.34). In particular, Lemma 3.6 follows 
immediately ■ 



Proof of Proposition 3.5. Start by multiplying (3.32) by e-h R(x a ,y a )da an( j then let s — > sq. Since 
P{x,y) = xyj^i j^iPi,j xl ~ 1 y : '~ 1 > see (1-4)) and since lim s ^. so x s y s eh R(^a,y a )da _ gee L emma 3 5 
we obtain that 

lim P(x s ,y s )efo R ( x ">y<*) da = 0, 
which concludes the proof of Proposition 3.5. ■ 
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Remark 5. When (xq, yo) G (0, l) 2 , we also have (x u , y u ) G (0, l) 2 for all u G (0, sq). Thus it is possible 
to plug approximations of the pij's and p^i's into (3.33) thanks to the terms (x u ) 1 and {y u V ■ Using 
that pi t i = 2d/(ri) + o(l/i) when i — > +oo, it is possible to find Iq sufficiently large so that: 



2 i>/ 7 ° 

/•So 

^ / V ((x u ) J + (y u y)e^ R (^v^du 
Jo i>/ 

= d / ( _5« + ) e / R(x a ,y a )da du 

J \l-x u l-y u J 



Thus: 

^0,2/0) = J 5>,i / i((x u y + (y u y)ef° R ^^ da du 
2 i=i J ° 

+ d f so ( x u (xl?-y u ) + Vuivi - vu) \f- R(Xa ,y a)dadu + o(x / 0+ i + y/o+i^ (336) 
Jo ^ ^ — x u 1 — y u J 

The latter expression shows that numerically, one can restrict to the computation of a finite number 
of probabilities p^i, for i ^ Iq- 

4 Numerical results 

In this section, we present two different ways of approximating the extinction probabilities pij. 
4.1 Probabilistic algorithm 

A first possibility, if we are interested in a given initial condition is to approximate Pij by 

Monte-Carlo simulations. For T > large, we simulate M paths (Xf , Y t )te{i,...,T} started at (i,j), 
for £ G {1,...,M}, independent and distributed as the process (^Q, Yt)t£{i,...,T}- The extinction 
probability is estimated by: 

1 M 

P M > T = mY1 %t<T,x|y/=o}- 
1=1 

The estimator pm,t is the proportion of paths that have gone extinct before time T. 



Proposition 4.1. Let be the initial condition. The estimator pm,t has the following properties: 

(i) It is a convergent and unbiased estimator of Pj j [tq ^ T] . 

(ii) Its variance is ¥ij[ro ^ T](l — Pjj[ro ^ T\)/M, and hence we have the following asymptotic 
95% confidence interval for pij : 



, ..... ,PM,t(1 ~Pm,t) ~ . p.,* Pm,t(1~Pm,t) 
p M ,T - 1.96a/ — - — —;pm,t + l-96i 



M ' V M 



(4.1) 



Proof. These results are straightforward consequences of the law of large numbers and central limit the- 
orem, given that lp is ^ x e Y £ =o} are independent Bernoulli random variables with parameter Pij[tq ^ 
T\. U 
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Computing the extinction probabilities by Monte-Carlo methods yields good results if we are 
interested in a given initial condition We then have a complexity of order M x T. However, 

biologists may be interested in investigating the extinction probabilities when the initial condition 
varies, and the method become computationally expensive. 



4.2 Deterministic algorithm 

For numerical approximations, we restrict ourselves to the computation of {pi,j)ije{l,...,N} f° r a positive 
(large) integer N. In this case, (1.7) can be approximated by the solution to a linear system. 

Let us define = (p 1; x, . . . ,Pi,n,P2,i, ■ ■ ■ ,P2,N, ■ ■ ■ ,PN,i, ■ ■ ■ ,Pn,n) T and T N is a iV 2 x 7V 2 -matrix 
with five non-zero diagonals: 



( 


A x 


D 





\ 






A 2 


D 
























D 


\ 










Bn,n-i An J 



where D = 2 (J+d) ^ JV ' ^ (* e {!>••■ > N}) and (i S {2, ... N}) are the N x A^-matrices 



A: 



d 2 
r+d i+2 



2(r+d) 
-1 

d 3 
r+d i+3 



B 



Li— 1 



d i 
r+d i+l 



d i 
r+d i+2 



d N 
r+d i+N 



\ 





r 

2(r+d) 
- 1 J 





d i 
r+d i+N 
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Let us also define the vector bjy = (6iat, . . . , &jvjv) T £ M, NxN such that: 



?ijv = - 



°NN 



d_ 

r+d 

_d 1_ 

r+d 1+2 



d 1 
r+d 1+(7V-1) 
d 1 _i_ r ~ 
\ r+dl+JV 2(r+d)^ 1 . Ar + 1 / 



d 1 
r+d 1+i 









farie{2,...,JV-l} J 



/ d 1 _i_ r ~ 
' r+dl+iV + 2^)™ + !-! 



2(^d)PN+l,N-l 
\ W+d) (PN+l,N +PN,N+l) J 



where Pi,jv+i an d Pn+i,n are approximations of Pi ; N+i given by Proposition 2.3. With these notations, 
(1.7) rewrites as 

TnPn = bN- 



4.3 Results 

We start with r = 3 and d = 2. For the Monte-Carlo simulation, we use M = 200 and T = 5000. For 
the deterministic method, we use ./V = 50, so that G {1, . . . , 50} 2 . Estimators of the extinction 
probabilities pf^ and pf^ obtained respectively from the methods of Sections 4.1 and 4.2 are plotted 
in Figure 5. The results given by both methods are very similar, as shown by the statistics of Table 1. 

pfj -pfj) 
For the latter, we consider only 



In Table 1, we compute the square difference between the two predictions (p\ L - —pf-) 2 -, the absolute 
difference \p^j 



p?)\ and the relative difference 



the couples where pf^ 



1,3 



whatever the value of pfj) 



and plj 



do not vanish (else, the fraction is either not defined or either 1 





Mean 




St.dev 


Min 




Max 


Square error 


3.24 10" 


-s 


2.47 10~ 4 


1.68 10" 


36 


4.63 10~ 3 


Absolute error 


1.34 10" 


-3 


5.53 10~ 3 


1.30 10- 


18 


6.81 10~ 2 


Relative error 


4.49 10" 


-2 


1.71 1Q- 1 


9.80 10" 


-3 


9.71 10- 1 



Table 1: Square, absolute and relative differences between the predictions of the stochastic method 
when r = 3 and d = 2 (Section 4-1) and of the deterministic methods (Section 4-. 2). Recall that with 
M = 200, the width of the confidence interval (4.1) is 6.92 10~ 2 . 



To carry further the comparison of the stochastic and deterministic method, and to observe the 
influence of N on the quality of the approximation, we compute the relative quadratic error 



^l<i,i<10 IPy Pi,j) y Z^l^iJ^lO [Pi,j Pi,j ) , . 

or : (4.2) 



{PiJ ) 
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(a) 



(b) 




Figure 5: Estimation of the extinction probabilities Pij's when r = 3 and d = 2: (a) with the 
Monte- Carlo method of Section 4-1- (b) with the deterministic method of Section 



when pi ■ is the deterministic approximation for N £ {10,..., 50} and pf- is either given by the 



deterministic approximation with N = 50, or by the stochastic approximation p^j with M = 200. 



p^j , the decrease in the quadratic errors stops around N = 18 around 



In the first case when p\ ■ 

0.0891. This corresponds roughly to the stochastic error of the law of large numbers (4.1) which 
depends only on M, In the second case, when p\ ■ is the deterministic approximation with TV = 50, 
the relative quadratic errors decrease exponentially fast in exp(-0.6842 N) (R 2 = 99.92%). 

In a second experiment, we choose r = 2.002 and d = 2. This case is more interesting in population 
ecology, since small populations are of interest when they are fragile and endangered species. For the 
Monte-Carlo simulation, we use M = 200 and T = 5000. For the deterministic method, we use 
N = 100, so that (i,j) £ {1, . . . , 100} 2 . The estimated extinction probabilities pj^ and p\j are 
plotted in Figure 6, and statistics are computed in Table 2. Again, results from both methods are 
similar. This is confirmed by computing the relative quadratic errors, with the p^J's obtained from 



d 3 ) - d 1 ) 



p\ '■ 's from the Monte-Carlo method. The decrease of this 



the deterministic method and the p, 
error is exponential with N in exp(— 0.0619 N) (R 2 = 98.98%) as shown in Figure 6(c). It can be 
noticed that in this case, the per! formances of the Monte-Carlo method match better the one of 
the deterministic algorithm. This is due to the fact that Monte-Carlo methods fail to produce good 
estimates of small probabilities (see [5] and references therein). 

When the probabilities p^J's are given by the deterministic method with iV = 50 in (4.2), we have as 
in the previous case (r = 3) an exponential decrease of the relative quadratic error in exp(— 0.1092 N) 
(R 2 = 95.07%). 





Mean 


St.dev 


Min 




Max 


Square error 


6.27 10~ 3 


6.25 10" 3 


8.42 10" 


10 


5.31 10" 2 


Absolute error 


6.89 10~ 2 


3.90 10- 2 


2.90 10" 


-5 


2.31 10" 1 


Relative error 


2.22 10" 1 


2.17 lO" 1 


1.02 10" 


-4 


1 



Table 2: Square, absolute and relative differences between the predictions of the stochastic method 
when r = 2.002 and d = 2 (Section 4-1) and of the deterministic methods (Section 4-%)- As in Table 
1, since M = 200, the width of the confidence interval (4.1) is 6.92 10 -2 . 
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(a) (b) 




Figure 6: Estimation of the extinction probabilities pij 's when r = 2.002 and d = 2: (a) with the 
Monte- Carlo method of Section 4-1- (b) with the deterministic method of Section 



Co' i viVitpnce oftho {p. .;|. ; .-. ,<■[(: for toryii!;; A" n-ith r '>. (102 ajuj d 2 




Figure 7: Estimation of the extinction probabilities pij 's when r = 2.002 and d = 2: Decrease with 
N of the log of the relative quadratic errors (4.2) between the deterministic and Monte-Carlo methods 

(m = 200;. 
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